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1.0  INTRODUCTION 

\ 

4 

Optical  signal  processing  makes  possible  rapid  coordinate  transformations, 
optical  pattern  recognition,  and  matrix-matrix  multiplication.  In  the  present 
contract,  DSI  has  demonstrated  several  significant  accomplishments.  Among  these 
are:  (a)  the  design  and  operation  of  a  spatial-frequency-multiplexed  coherent 
optical  processor;  (b)  the  application  of  outer  product  operation  to  matrix- 
matrix  multiplication;  and  (c)  the  description  of  three  methods  for  performing 
Hankel  transforms  with  optical  or  digital  processors. 
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2.0  RESULTS 


The  accomplishments  described  below  are  important  in  a  number  of  areas. 

The  optical  information  processing  projects  described  have  application  to  mode 
analysis  of  high  energy  laser  resonators,  alignment  of  space  telescopes,  and 
tomographic  reconstruction.  The  progress  in  coherent  optical  moment  generation 
is  of  value  in  image  data  compression  and  optical  pattern  recognition  of 
statistical  images.  Finally,  the  acoustooptic  matrix  signal  processor  con¬ 
structed  is  useful  in  radar  signal  processing,  passive  sonar  surveillance,  and 
adaptive  beam  forming. 

2.1  MULTIPLEXED  COHERENT  OPTICAL  PROCESSOR 

A  coherent  optical  processor  designed  for  calculating  generalized  moments 
of  a  two-dimensional  pattern*  is  described  in  Enclosure  (1).  In  order  to  provide 
for  parallel  computation  of  multiple  moments,  a  spatial-frequency  multiplexing 
scheme  was  used.  A  computer-generated  holographic  mask  was  shown  to  provide 
complete  flexibility  in  choosing  moment  generated  functions.  The  calculation 
of  five  geometric  moments  was  experimentally  demonstrated  for  simple  objects. 
These  geometric  moments  correspond  to  x,  y,  xy,  x^,  and  y2. 

A  computer-generated  holographic  mask  was  constructed  to  compute  the  geo¬ 
metric  moments  corresponding  to  the  generating  functions,  x,  y,  xy,  x^,  and  y^. 

In  order  to  reduce  the  dynamic  range  problems,  when  the  hologram  was  construct¬ 
ed  the  scale  of  the  generating  functions  was  chosen  such  that  the  linear  and 
quadratic-moment  values  would  be  equal  for  a  symmetric  object  threefourths  the 
linear  size  of  the  mask.  The  spacing  chosen  avoided  crosstalk  between  adjacent 
moments.  Since  several  generating  functions  are  encoded  on  a  single  mask,  the 
optical  processor  is  able  to  calculate  several  moments  in  parallel. 

2.2  APPLICATION  OF  OUTER  PRODUCT  OPERATION  TO  MATRIX-MATRIX  MULTIPLICATION 

The  systolic  architecture  proposed  by  Caulfield  and  Rhodes^  for  matrix- 
matrix  multiplication  was  utilized.  This  approach  bypasses  the  two-dimensional 
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real  time  spatial  light  modulators  required  by  other  architectures.  Various 
implementations  of  this  optical  processor^  can  use  one  of  the  following  modula¬ 
tion  mechanisms:  (a)  electrooptic,  (b)  direct  driven  LED  array,  and  (c)  acousto¬ 
optic  (e.g.,  Bragg  cells).  These  three  approaches  are  described  in  Enclosure  (2). 

2.3  OPTICAL  AND  DIGITAL  TECHNIQUES  FOR  PERFORMING  HANKEL  TRANSFORMS 

Generalized  Hankel  transforms  are  useful  in  analyzing  the  effect  of  cir¬ 
cularly  symmetric  optical  systems  on  arbitrary  inputs.  Some  examples  of  such 
systems  are  complex  laser  resonators  and  space  telescopes.  Three  methods  for 
performing  Hankel  transforms  with  optical  or  digital  processors  were  described^. 
The  first  method  was  applicable  when  the  input  data  is  available  in  cartesian 
(x-y)  format  and  used  the  close  connection  between  generalized  Hankel  transform 
and  the  two-dimensional  Fourier  transform  in  cartesian  coordinates.  The  second 
method  was  useful  when  the  input  data  is  in  polar  (r-0)  format  and  used  change 
of  variables  to  perform  the  nc^  order  Hankel  transform  as  a  correlation  inte¬ 
gral.  The  third  method  utilized  the  von  Neumann  addition  theorem  for  Bessel 
Functions  to  extract  the  Hankel  coefficients  from  a  correlation  between  the 
radial  part  of  the  input  and  a  Bessel  function.  Initial  experimental  results^ 
obtained  for  optical  implementation  of  the  first  two  methods  are  presented 
in  Enclosure  (3). 

The  analysis  of  complex  optical  systems  is  greatly  facilitated  by  two- 
dimensional  Fourier  techniques.  The  effect  of  an  optical  system  on  arbitrary 
inputs  is  easily  described  by  a  transfer  function  in  the  Fourier  domain.  Gen¬ 
eralized  Hankel  transform  is  similarly  useful  when  dealing  with  a  circularly 
symmetric  (or  axisymmetric)  system  for  arbitrary  inputs.  This  situation  is 
encountered  in  performing  mode  analysis  on  the  output  of  a  slightly  misaligned 
laser  resonator  as  well  as  in  aligning  space  telescopes.  An  optical  method  for 
performing  mode  analysis  via  generalized  Hankel  transform  has  the  unique  advan¬ 
tage  of  preserving  the  phase  of  the  wavefront  to  be  analyzed. 
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A  coherent  optical  processor  capable  of  calculating  generalized  momenta  of  a  two-dimensional  pattern  is  described. 

A  spatial-frequency  multiplexing  scheme  is  used  to  provide  for  parallel  computation  of  multiple  moments.  The 
use  of  a  computer-generated  holographic  mask  permits  complete  flexibility  in  choosing  moment-generating  func¬ 
tions;  e.g.,  the  functions  could  be  complex  or  have  a  predetermined  weighting  function.  Experimentally,  calcula¬ 
tion  of  five  geometric  moments  (corresponding  tox,y,xy,x2,  and  y})  is  demonstrated  for  simple  objects.  The  spe¬ 
cial  features  of  the  proposed  coherent  optical  processor  and  its  space-bandwidth  requirements  are  also  discussed. 


The  role  of  moments  and  moment  invariants  as  global 
features  in  pattern  recognition  is  well  known.1-3  Co¬ 
herent  optical  processors,  for  calculating  geometric 
moments  of  images  have  been  recently  proposed,4'5  but 
both  of  these  processors  are  restricted  to  calculating 
geometric  moments  only.  In  addition,  the  processor  of 
Ref.  4  is  rather  complicated,  involving  interferometric 
setups;  the  processor  of  Ref.  5  has  an  added  factor  1/ 
ip'.ql)  associated  with  the  pq th  moment,  further  com¬ 
pounding  the  dynamic-range  problems.  In  t  his  Letter 
we  propose  a  spatial-frequency-multiplexed  coherent 
optical  processor  that  overcomes  these  limitations.  We 
describe  the  design  and  operation  of  this  processor, 
present  the  initial  experimental  results  obtained  with 
it,  and  summarize  its  special  features  and  limitations. 

The  generalized  moment  Cmn  of  a  two-dimensional 
pattern  fix,  y)  is  defined  as 

Gmn  58  x:x:  f(x,y)g  mn  (x,y)dxdy,  (1) 

where  gmn (*,  y)  is  the  generating  function.  One  way 
of  calculating  the  inner  product  of  Eq.  (1)  optically  is 
to  multiply  the  input  fix,  y )  by  the  generating  function 
gmnix,  y)  and  then  take  a  two-dimensional  Fourier 
transform  of  the  product  The  moment  Gm„  is  obtained 
at  the  origin  of  the  frequency  plane: 

Gmn  *  x:x:  fix,y)g  mn  (f  i  3') 

X  exp(— y(ux  +  oy)]dxd.v  u_0  .  (2) 

1—0 

If  gmn  (x,  y)  is  multiplied  by  a  spatial  carrier  of  the  form 
exp[/(muoX  +  nt/oy)],  then  the  desired  moment  is  ob¬ 
tained  at  u  *  muo,  v  *  nvo  in  the  frequency  plane 

Gmn  *  x:x:  /U,y)gmn(x,y)exp)-;'((i/  -  mu{ ,)x 


Now  if  the  function  gix,  y ),  which  multiplies  the  input 
fix,  y),  has  the  form 
M  N 

gix,y)  m  £  L  gnnix,y)exp\jimuox  +  m>oy)], 

m  M  n--JV 

(4) 

then  the  output  of  the  optical  processor  will  contain  the 
desired  generalized  moments  Gm„  (m  *  -Af, ...  ,  0, . . . , 

M  and  n  =  -N, . . . ,  0 . N)  provided  that  the  unit 

carrier  frequency  (u0,  Do)  satisfies  the  conditions  dis¬ 
cussed  below.  Figure  1  is  a  schematic  diagram  of  the 
optical  processor.  It  consists  simply  of  a  mask  con¬ 
taining  gix,  v)  described  by  Eq.  (4)  and  a  Fourier- 
transform  lens.  A  two-dimensional  array  of  detectors 
measures  the  light  intensity  at  discrete  points  in  the 
Fourier  plane,  giving  |  Gmn ( 2-  Since  the  function  gix, 
y)  is  bipolar  in  the  case  of  geometric  moments  and  could 
be  complex  for  some  other  choices  of  generating  func¬ 
tions,  it  is  necessary  to  make  the  mask  holographically. 
In  particular,  a  computer-generated  hologram  will 
provide  maximum  operational  flexibility  to  the  coherent 
optical  processor. 

To  demonstrate  the  operation  of  the  coherent  optical 
processor  described  above,  a  computer -generated  ho¬ 
lographic  mask  was  constructed  to  compute  the  geo¬ 
metric  moments  mm,  moi,  mu,  m20,  and  mo2  (corre¬ 
sponding  to  generating  functions  x,  y,  xy,  x2,  and  y2, 
respectively).  When  the  hologram  was  constructed  t.he 
generating  functions  were  scaled  such  that  the  linear- 
and  quadratic-moment  values  would  be  equal  for  a 
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Fig.  2.  Plot  of  the  generating  functions  x  and  x'1  with  x 
multiplied  by  a  weight  of  0.5. 


Fig.  3.  Arrangement  of  the  five  geometric  moments  in  the 
output  plane. 

symmetric  object  three  fourths  the  linear  size  of  the 
mask  (Fig.  2).  This  was  done  to  alleviate  the  dy- 
namic-range  problems  in  making  the  mask.  The  ar¬ 
rangement  of  the  five  moments  in  the  output  plane  is 
shown  in  Fig.  3.  The  spacing  was  chosen  to  avoid  cross 
talk  between  adjacent  moments.  It  should  be  noted 
that  the  arrangement  is  determined  by  the  spntinl  car¬ 
rier  multiplying  the  different  generating  functions  and 
hence  is  completely  under  experiment  n I  emit  ml.  The 
computer-generated  holograms,  which  used  the  Lee- 
Burckhardt  encoding  scheme,8  contained  216  X  210 
pixels  and  had  physical  dimensions  of  10.2  nun  X  10.2 
mm.  The  output  was  detected  by  a  video  camera  with 
a  PbO  target.  The  y  of  the  camera  was  adjusted  to  be 
»1;  hence  the  electrical  signal  was  proportional  to  the 
light  intensity.  The  output  of  the  camera  was  digitized 
and  stored  in  a  minicomputer  for  further  analysis. 

The  coherent  optical  processor  was  used  to  calculate 
these  five  geometric  moments  for  simple  inputs.  A 
photograph  of  the  output  for  a  binary  square  input 
image  is  displayed  as  Fig.  4.  The  moments  are  mea¬ 
sured  at  the  center  of  each  of  t  he  five  pat  t  erns.  The 
and  moi  moments  are  seen  to  lie  zero,  indicating  that  the 
input  had  symmetry  about  the  V  and  .V  axes,  respec¬ 
tively.  A  line  scan  through  the  ni„,  pattern  clearly 


shows  the  null  at  the  center,  as  depicted  in  Fig.  5(a). 
When  the  input  was  translated  by  ~10%  of  its  lineur 
dimension  along  the  X  axis  (thus  removing  the  sym¬ 
metry  about  the  Y  axis),  a  nonzero  value  for  mm  is  ob¬ 
tained,  as  is  seen  in  Fig.  5(b). 

The  line  scans  through  m2o  and  m<)2  moments  of  a 
binary  rectangular  input  with  a  2:1  aspect  ratio  and  X- 
and  V-axis  symmetry  are  displayed  in  Fig.  6.  The 
theoretical  value  of  the  ratio  |m2o/mo2|2  is  16,  whereas 
the  experimentally  measured  value  is  16.3,  corre¬ 
sponding  to  an  error  of  ~1.9%.  A  detector  with  better 
dynamic  range  and  linearity  should  result  in  improved 
accuracy. 

In  this  Letter  we  have  described  a  coherent  optical 
processor  designed  for  calculating  generalized  moments 
of  a  two-dimensional  pattern.  A  computer-generated 
hologram  is  used  to  encode  several  generating  functions 
on  one  mask  with  spatial-frequency  encoding.  This 
special  feature  enables  the  processor  to  calculate  several 
moments  in  parallel.  The  use  of  a  computer-generated 
hologram  makes  it  possible  to  encode  complex  gener¬ 
ating  functions  and  permits  predetermined  weighting 
functions  to  be  used  for  dynamic-range  consider¬ 
ation. 


Fig.  4.  Photograph  of  the  output  of  the  optical  processor. 
The  input  was  a  square  binary  object  with  x  and  y  sym¬ 
metry. 


Fig.  5.  (n>  Line  scan  through  the  m  m  moment  of  a  binary 
square  with  x  symmetry,  (b)  Line  scan  through  the  mm 
moment  of  a  binary  square  shifted  in  the  x  direction  wit  h  re¬ 
spect  to  its  position  in  (a). 
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(a)  (b) 

» 

Fig.  6.  Line  scans  through  two  geometric  moments  of  a 
rectangle  with  2:1  aspect  ratio:  (a)  mw,  (1>)  m»<. 


IB,  *  Bml 

Fig.  7.  Distribution  in  the  output  plane  when  i/,i  =  H,  + 
leading  to  zero  cross  talk  between  adjacent  momenls. 


The  maximum  number  of  moments  that  can  be  cal¬ 
culated  in  parallel  is  determined  by  the  single-sided 
bandwidth  of  the  system  (including  the  computer¬ 
generated  hologram)  Rs,  the  single-sided  bandwidth  of 
the  input  funct  ion  B,,  and  the  single-sided  bandwidth 
of  the  generating  function  Bm.  Since  the  output  con¬ 
tains  a  convolution  between  the  Fourier  transform  of 
the  input  and  the  Fourier  transform  of  the  generating 
function,  each  term  will  be  extended  over  a  region  2{B, 
+  Ii„, )  in  the  frequency  plane.  The  spacing  of  the  mo¬ 
ments  (determined  by  u o  in  u  direction  and  by  uo  in  u 
direction)  should  be  such  that  at  points  (mito,  m>o)  only 
the  Gmn  will  be  nonzero;  i.e.,  uo  =  t>o  =  (5/  +  Bm),  as  is 
shown  in  Fig.  7  for  one  dimension.  Therefore  the  total 
number  of  moments  that  can  be  calculated  in  parallel 
is  given  by 


N 


o 


'  Bs  \ 
Bj  +  Bm) 


(5) 


This  number  can  be  increased  substantially  by  using 
spectrum-shaping  techniques  on  the  input;  these 
techniques  will  be  discussed  in  a  subsequent  publica¬ 
tion. 
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The  use  of  generating  functions  other  than  xn'y" 
could  offer  significant  advantages.  For  given  classes 
of  functions,  certain  combinations  of  geometric  mo¬ 
ments  can  be  determined  to  be  optimum  for  minimizing 
intraclass  spread  while  maximizing  interclass  separation 
by  using  standard  algorithms."  These  combinations 
can  then  be  encoded  directly  in  the  computer-generated 
holographic  mask,  thus  enabling  the  coherent  optical 
processor  to  carry  out  classification  of  statistical  pat¬ 
terns.  To  extract  complex-value  moments  in  the  out¬ 
put,  the  four  quadrature  components  could  be  treated 
independently  and  measured  ns  positive  real  num¬ 
bers. 
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Optical  processors  for  multiplying  two  matrices  have  been 
described  in  the  literature.1-4  Considerable  work  has  been 
reported  recently  on  performing  vector-matrix  multiplica¬ 
tions  with  incoherent  optical  processors.5-15  It  is  suggested 
that  matrix-matrix  multiplications  can  be  implemented  with 
these  processors  by  inputting  one  matrix  a  row  at  a  time  while 
the  second  matrix  is  encoded  on  a  2-D  mask.  The  output 
matrix  is  then  available  one  row  at  a  time.  All  these  schemes 
have  one  drawback  in  that  they  require  at  least  one  of  the 
input  matrices  on  a  2-D  mask,  thus  making  it  difficult  to 
update  that  matrix  at  high  frame  rates  with  currently  available 
2-D  spatial  light  modulators  (SLM).  An  optical  processor 
based  on  a  systolic  architecture  was  recently  proposed  To 
overcome  this  drawback."  But  such  a  processor  requires  2N 
—  1  parallel  channels  to  handle  N  X  N  matrices  and  tnkes  4N 
clock  cycles  to  perform  one  vector-matrix  multiplication. 
Thus  it  does  not  fully  utilize  the  parallel  processing  capabil¬ 
ities  of  optical  systems. 

These  processors  are  all  based  on  performing  the  inner 
product  between  two  vectors  for  generating  one  element  of  the 
output  vector  (or  matrix).  Thus  they  use  dimensionality- 
reducing  operations  (from  N  to  1).  An  alternate  approach 
would  involve  dimensionality-increasing  operations  such  as 
an  outer  product  between  two  vectors  (from  N  to  N-).  Pro¬ 
posals  have  been  made  to  use  the  outer  product  between  two 
vectors  to  generate  the  covariance  matrix  of  those  vectors  8  8 
In  this  Letter  we  describe  the  application  of  the  outer  product 
operations  to  matrix-matrix  multiplication. 

The  matrix  multiplication  between  two  N  X  N  matrices  A 
and  B  can  be  stated  as  follows10: 

C  =  A  B  (1) 


LIGHT  SOURCE  CROSSEO  PAIR  OF  TWO-DIMENSIONAL 

LINEAR  EO  MODULATORS  DETECTOR  ARRAY 


COLUMNS  ROWS  OUTER  PRODUCT 

OF  A  OF  £  MATRIX 


Fig.  1 .  Schematic  diagram  of  an  KO-EO  processor  for  matrix  mul¬ 
tiplication.  The  crossed  pair  of  linear  modulators  and  detector  plane 
are  shown  separated  for  the  sake  of  clarity. 
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Fig.  2.  Schematic  diagram  of  a  LED-EO  processor  for  matrix  mul 
tiplication.  Imaging  and  focusing  optics  are  omitted. 


Thus  the  ijth  element  of  C  is  given  by  the  inner  product  be¬ 
tween  the  ith  row  vector  of  A  and  the yth  column  vector  of  11. 
Alternatively  the  output  matrix  C  can  he  expressed  ns  a  sum 
of  N  matrices: 

C  *  L  C-  (2) 

■-1 

where 


Here  each  matrix  term  C’  in  the  summation  is  seen  to  corre¬ 
spond  to  an  outer  product  between  the  ith  column  vector  of 
A  and  ith  row  vector  of  H.  The  ooter  product  between  two 
vectors  can  he  performed  optically  by  crossing  two  linear-array 
light  modulators  as  described  in  Ref.  8.  The  summation  of 
the  different,  matrix  terms  in  Kq.  (2)  ran  be  performed  in  a  2-D 
integrating  detector  array.  To  perform  the  operation  of 
matrix  multiplication  between  two  N  X  N  matrices  A  and  B, 
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the  columns  of  A  and  rows  of  B  tire  applied  sequentially  to  the 
appropriate  elements  of  the  crossed  pair  of  linear  arrays  of 
light  modulators.  The  light  transmi  tted  through  the  crossed 
light  modulators  contains  the  outer  product  between  two 
vectors  which  is  detected  and  summed  in  a  2-1 )  detector  array. 
Thus  it  takes  N  clock  cycles  to  carry  out  the  matrix  multipli¬ 
cation  which  involves  N3  multiplies  and  A/:l  addil  ions.  The 
optical  processor,  therefore,  performs  /V*  multiplies  and  N- 
additions  in  parallel.  In  this  algorithm  both  of  the  input 
matrices  are  required  one  column  (row)  vector  at  a  time,  and 
the  well -developed  technology  of  1-D  light  modulators  can  he 
utilized  in  building  the  desired  optical  processor. 


the  acoustic  wave  in  the  Bragg  cell  in  time.  The  columns  of 
A  and  rows  of  B  modulate  an  rf  carrier  which  drives  the  ap¬ 
propriate  Bragg  cells.  When  all  the  elements  of  the  column 
(row)  vector  are  loaded  in  the  Bragg  cell,  the  light  source  is 
pulsed  and  the  outer  product  is  stored  in  the  detector 
array. 


The  choice  of  n  particular  implementation  of  the  processor 
will  he  guided  by  the  format  of  input  matrices  as  well  as  the 
constraints  imposed  by  the  volume  and  power  consumption 
requirements  of  specific  applications. 


7 

The  various  implementations  of  this  optical  processor  can 
use  one  of  the  following  modulation  mechanisms  for  the  linear 
array  of  light  modulators:  electrooptic  (EO)  modulations, 
direct  driven  LED  array,  or,  acoustooptic  (AO)  Bragg  cells. 


In  this  Letter  we  describe  three  implementations: 

(1)  Encoding  of  both  A  and  B  via  EO  modulators  (EO-EO 
processor)  leads  to  a  very  compact  optical  system  which  does 
not  require  any  imaging  or  focusing  components  such  as  lenses 
or  fibers.  The  schematic  diagram  of  the  optical  processor  is 
shown  in  Fig.  I,  where  the  two  crossed  light  modulator  arrays 
and  detector  are  shown  separated  for  the  sake  of  clarity.  The 
EO  modulators  can  be  constructed  out  of  PLZT  elect  rooptic 
ceramics  which  possess  good  switching  speed  and  require 
moderate  voltages.11 

(2)  The  schematic  diagram  of  the  optical  processor  where 
A  is  encoded  via  direct  modulation  of  a  LED  array  and  B  is 
encoded  via  EO  modulators  (LED-EO  processor)  is  shown  in 
Fig.  2.  Here  an  optical  system  spreads  the  LED  array  output 
in  the  horizontal  direction  while  focusing  in  the  vertical  di¬ 
rection.  The  linear  array  of  EO  modulators  can  consist  of 
point  modulators  instead  of  the  elongated  finger  modulators 
required  by  the  EO-EO  scheme.  The  emerging  light  is  then 
imaged  in  both  horizontal  and  vertical  directions  on  a  2-D 
detector  array 

(3)  The  optical  processor  based  on  encoding  bot  h  A  and  B 
via  Bragg  cells  (AO-AO  processor)  is  similar  to  the  crossed 
Bragg  cell  processors  described  in  the  literature  for  generating 
Woodward’s  ambiguity  functions12  (see  Fig.  3).  The  differ¬ 
ence  in  this  case  is  that  the  light  source  is  pulsed  so  ns  to  freeze 


PULSED  LIGHT 
SOURCE 


CROSSED  PAIR  OP 
BRAGG  CELLS 


TWO-DIMENSIONAL 
DETECTOR  ARRAY 


COLUMNS  OP  ROWS  OP  OUTER  PRODUCT 

A  B  MATRIX  C. 


Fig  3.  Schematic  diagram  of  an  AO-AO  processor  for  matrix  mul¬ 
tiplication.  Imaging  optics  as  well  as  spatial  filtering  are  omitted. 


The  overall  performance  of  this  optical  processor  will  be 
determined  by  the  operating  characteristics  of  the  light  source, 
the  light  modulators,  and  the  2-D  detector  array.  Since  the 
technologies  pertaining  to  the  light  sources  and  I-D  modu¬ 
lators  are  well-developed,  the  limitation  on  the  accuracy  and 
throughput  rate  of  the  processor  will  be  imposed  by  the  2-D 
detector  array.  The  main  choices  for  the  detector  array  are 
going  to  be  CCD  array,  photodiode  array,  and  vidicon  TV 
camera.  As  a  representative  example,  if  one  considers  the 
third  implementation  of  the  optical  processors  (AO-AO)  with 
the  vidicon  TV  camera  as  the  detector,  it  is  possible  to  mul¬ 
tiply  two  positive  1000  X  1000  matrices  in  30  n;  ec  corre¬ 
sponding  to  3  X  1010  multiplies  and  3  X  10l°  adds/sec  with 
8-bit  accuracy.  Such  a  matrix-matrix  multiplier  would  be 
useful  in  2-D  mathematical  transforms,  matrix  inversion 
problems,  and  pattern  recognition  among  other  tasks. 


Stimulating  discussions  with  D.  Stillwell,  J.  N.  Lee,  and  A. 
D.  Fisher  arc  grntefully  acknowledged. 
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Integral  transforms  involving  Bessel  function  kernels  are  useful  in  analyzing  effects  of  circularly  symmetric  optical 
systems  on  arbitrary  inputs.  Methods  for  performing  the  integral  transforms  optically  are  divided  into  two  catego¬ 
ries.  The  first  category  involves  input  data  available  in  cartesian  U,  y)  format  and  uses  the  close  connection  be¬ 
tween  the  desired  integral  transform  and  the  two-dimensional  Fourier  transform  in  cartesian  coordinates.  The 
second  category  involves  input  data  in  polar  (r,  0)  format  and  uses  methods  such  as  change  of  variables  to  perform 
the  integral  transform  as  a  correlation  integral.  Experimental  results  obtained  with  optical  implementation  for 
these  two  categories  are  presented. 


The  analysis  of  complex  optical  systems  is  greatly  fa¬ 
cilitated  by  two-dimensional. Fourier  transform  tech¬ 
niques.  The  effect  of  an  optical  system  on  arbitrary 
inputs  is  easily  described  by  a  transfer  function  in  the 
Fourier  domain.  Integral  transforms  involving  Bessel 
functions  are  similarly  useful  when  dealing  with  circu¬ 
larly  symmetric  (or  axisymmetric)  systems  for  arbitrary 
inputs.1  This  problem  is  encountered  in  the  mode 
analysis  on  the  output  of  a  slightly  misaligned  laser 
resonator23  or  in  the  alignment  of  space  telescopes.  An 
optical  method  for  performing  mode  analysis  via  inte¬ 
gral  transforms  will  have  the  unique  advantage  of  pre¬ 
serving  the  phase  of  the  wave  front  to  be  analyzed,  in 
contrast  to  the  more  usual  methods  involving  mea¬ 
surement  of  the  intensity  of  the  wave  front. 

It  is  well  known  that  when  a  two-dimensional  func¬ 
tion  has  circular  symmetry,  its  Fourier  transform  is  also 
circularly  symmetric.  It  can  be  shown  that  in  such  a 
case  the  Fourier  transform  is  equivalent  to  the  Hankel 
transform  of  the  input,  which  is  defined  as 

F°(p)  =  f{r)Jo(pr)rdr,  (1) 

where  Jo(pr )  is  the  zeroth-order  Bessel  function. 1  Thus 
in  dealing  with  circularly  symmetric  systems,  the 
Hankel  transform  (which  is  a  one-dimensional  integral) 
can  be  used  instead  of  the  two-dimensional  Fourier 
transform  if  the  inputs  are  aiso  circularly  symmetric. 
This  technique  can  be  extended  to  functions  of  the  form 
f(r )  expijnd)  (n-fold  symmetric)  by  using  the  Hankel 
transform  of  order  n  defined  as 

Fn(p)  =  f(r)Jn(pr)rdr,  (2a) 

where  Jn  ( pr )  is  the  nth-order  Bessel  function.  In  this 
case  the  two-dimensional  Fourier  transform  with  re¬ 
spect  to  (x,  y) — indicated  by  FTo — is  related  to  /*’„(/>) 
by1 

FT^lfir)  exp(/n0)]  =  2ir F„(p)  exp(-y>i</>),  (2b) 
where  ( p ,  <t> )  are  polar  variables  in  the  Fourier  plane. 


The  generalized  Hankel  transform,  Fnn(p),  can  be 
defined  for  an  arbitrary  function  f(r,  d)  as  follows1: 

f(r,  6)  -  f;  fn(r)  exp(jnd),  (3a) 

nm— <* 

Fnn(p)  =  /n(r)t/„(pr)rdr.  (3b) 

The  generalized  Hankel  transform  thus  involves  a 
Fourier  series  expansion  of  f(r ,  6)  with  respect  to  d, 
followed  by  the  operation  of  the  Hankel  transform  of 
order  n  on  fn(r).  In  the  following  we  will  describe 
techniques  and  experimental  results  for  optical  imple 
mentation  of  the  generalized  Hankel  transforms. 
These  various  techniques  are  each  applicable  in  dif¬ 
ferent  circumstances. 

In  optical  processors  a  two-dimensional  Fourier 
transform  with  respect  to  the  cartesian  coordinates  (x, 
v)  is  performed  very  easily  with  the  help  of  a  simple 
spherical  lens.4  The  equivalence  of  the  two-dimen¬ 
sional  Fourier  transform  and  the  Hankel  transform  for 
circularly  symmetric  functions  was  noted  above.  For 
arbitrary  functions,  f(x,  y),  the  following  relation  ex¬ 
ists: 

FT>[f(x,  y)]  =  (7(p,  <t>)  -  2ir  £  Fnn(p)  exp(~jn<p), 

n  = — 

(4) 

where  7(p,  <t>)  is  the  two-dimensional  Fourier  transform 
of  the  input  in  terms  of  the  polar  coordinates  (p,  <t>). 
The  generalized  Hankel  coefficients  Fnn(p)  can  thus  be 
extracted  by  performing  a  one-dimensional  Fourier 
transform  on  3{p,  tf>)  with  respect  to  <t>, 

F„n(p)  =  J*  3(p,  <fi)  exp(jnd>)dct>.  (5) 

Therefore,  a  spherical  lens  first  performs  the  two- 
dimensional  Fourier  transform  with  respect  to  the 
cartesian  coordinates  on  the  input.  A  suitably  designed 
computer-generated  hologram  then  performs  the 
coordinate  transformation  (cartesian  to  polar)  on  the 
Fourier  transform,  generating  7(p,  <t>).h  This  is  followed 
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Fig.  1.  Experimental  results  for  zeroth-order  Hnnkel 
transform,  Fo(p),  of  Hr  -  a),  (a)  Photograph  of  the  output, 
(b)  A  plot  of  |-7o(ap)|2,  which  is  the  theoretically  expected 
result,  (c)  A  line  scan  through  the  origin  of  the  pattern  in  (a), 
giving  |/'o(p)|-  versus  p. 


Fig.  2.  Experimental  results  for  first-order  Hankel  trans¬ 
form,  F\{p),  of  5  (r  -  a),  (a)  A  plot  of  |Ji(ap)|-  versus  p.  which 
is  the  theoretically  expected  result,  (h)  A  line  scan  through 
the  origin  of  the  two-dimensional  Fourier  transform,  giving 
)  Fi(p)  | 2  versus  p. 

by  a  cylindrical  lens  to  implement  the  one-dimensional 
Fourier  transform  with  respect  to  <t>. 

Initial  optical  experiments  established  the  connection 
between  a  two-dimensional  Fourier  transform  with 
respect  to  cartesian  coordinates  and  a  Hnnkel  t  ransform 
of  order  n  for  a  function  of  the  form  fir)  exp  (JnO)  (Eq. 
(2b)].  The  optical  system  used  was  a  standard  two- 
dimensional  Fourier-transform  arrangement.'  A 
computer-generated  hologram  was  used  to  encode  a 


expO'nfl)  dppendence  onto  the  input.  The  output  was 
detected  by  a  TV  camera,  which  measures  the  light  in¬ 
tensity  in  the  Fourier  plane  of  the  input.  The  input 
used  in  these  experiments  had  an  r  dependence  given 
by  6(r  -  a),  thus  corresponding  to  a  thin  ring  of  radius 
a  in  the  cartesian  (x,y)  plane,  The  nth-order  Hankel 
transform  of  5(r  -  a)  is  a  Jn(ap),  thus  allowing  easy 
comparison  with  experiments.  Figure  1  shows  the  re¬ 
sults  corresponding  to  the  Hankel  transform  (i.e.,  n  = 
0  and  exp(/n0)  =  1].  In  Fig.  1(b)  the  function  \Jo{ap)\2 
is  plotted,  which  is  then  compared  with  a  line  scan 
through  the  origin  |Fig.  1(c)]  of  the  output  shown  in  Fig. 
1(a),  i.e.,  <t>  -  0.  Very  good  qualitative  agreement  be¬ 
tween  the  theoretical  and  experimental  results  is  ob¬ 
tained.  Figures  2(a)  and  2(b)  present  the  results  for  the 
Hankel  transform  of  order  one  for  the  same/(r),  i.e.,  Hr 
-  a )  [but  here  the  9  dependence  is  exp(jO)].  Again  good 
qualitative  agreement  was  seen  between  theory  and 
experiment.  This  system  also  performed  the  Hankel 
transform  of  order  two  on  the  same  input  with  good 
results.6 

If  the  input  is  polar  ( r,d )  formatted,  the  direct  ap¬ 
proach  given  in  Eq.  (3)  has  to  be  followed  in  obtaining 
the  generalized  Hankel  transform.  The  first  part  of  the 
operation,  which  involves  a  Fourier  series  expansion  in 
variable  0,  requires  the  r  and  9  coordinates  be  mapped 
along  orthogonal  axes.  Then  Eq.  3(a)  is  performed 
optically  using  a  cylindrical  lens.  The  calculation  of  the 
Hankel  transform  of  order  n  of  /„(r),  the  nth  coefficient 
of  expansion,  is  less  straightforward,  since  it  involves 
optically  performing  a  space-variant  operation;  hence, 
this  was  the  main  aim  here. 

One  way  of  converting  a  space-variant  operation  into 
a  shift-invariant  operation  is  to  employ  appropriate 
change  of  variables.  In  the  case  of  the  Hankel  trans¬ 
forms  of  order  n,  the  following  procedure  was  described 
by  Siegman  for  implementing  the  space-variant  oper¬ 
ation  as  a  correlation  integral  on  a  digital  processor.7 
From  the  definition  (2a)  one  obtains 

Fn{y)  =  §  f(x)Jn(x  +  y)  dx,  (6) 

where r  =  r0exp(ax),p  =  p0exp(ay),Pa(y)  *  pFn(p), 
t(x)  =  rf{r),  3n(x  +  y)  =  arpjn(rp).  The  algorithm, 
therefore,  consists  of  first  linearly  weighting  the  input 
f(r)  and  performing  the  r  -*  x  coordinates  transform. 
This  new  input  is  then  correlated  with  a  similarly 
weighted  and  coordinate  transformed  nth-order  Bessel 
function  to  give  the  desired  Hankel  transform  also  in 
linearly  weighted  and  coordinate-distorted  form.  In 
any  physical  system  the  correlation  integral  will  be 
performed  over  a  finite  interval,  giving  rise  to  truncation 
errors.  Also,  if  the  input  is  sampled  in  the  x  domain, 
the  sampling  rate  should  be  adequate  to  represent  the 
function  accurately  in  the  x  domain.7  8 

Since  the  operations  of  coordinate  transformation 
and  correlation  can  be  performed  by  an  optical  pro¬ 
cessor,  an  optical  system  can  be  designed  and  used  to 
calculate  Fnn(P )•  Computer-generated  holograms  can 
be  used  to  perform  r  —  x  coordinate  transformation  as 
well  ns  t  o  encode  the  Fourier  plane  filter  with  impulse 
response  Jn  (x )  in  a  one-dimensional  correlator.  The 
second  dimension  of  the  optical  processor  can  be  used 
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to  perform  different-order  Hankel  transforms  on  dif¬ 
ferent  inputs,  thus  achieving  multichannel  operation. 

In  the  experiments,  both  the  input  and  the  Fourier 
plane  filter  were  encoded  by  computer-generated  ho¬ 
lograms.  The  linear  weighting  and  the  coordinate 
transformation  of  the  input  was  performed  by  the  dig¬ 
ital  computer  before  doing  the  holographic  encoding. 
The  computer-generated  holograms  used  the  Lee- 
Burckhardt  technique9  and  contained  128  pixels.  The 
optical  system  is  depicted  in  Fig.  3.  This  system  was 
•then  used  to  perform  the  Hankel  transform  on  two 
different  inputs /i(r)  =  p\J0(p\r)  and  /o(r)  =  p  *Ju(p>r). 
The  Hankel  transform  of  fi(r)  and  f2(r)  should  be  f>(p 
—  p\)  and  b(p  -  pv),  respectively.  The  results  of  a 
computer  simulation  of  this  algorithm  are  depicted  in 
Fig.  4.  The  finite  width  of  the  peak  and  the  appearance 
of  sidelobes  are  due  to  the  finite  limits  of  integration. 
To  facilitate  easy  comparison  with  the  experimental 
results  IFoCy)!2  was  plotted  versus  y  instead  of  F()(/>) 
versus  p.  Figure  5  shows  the  experimental  results  ob¬ 
tained.  The  optical  output  was  detected  by  a  1024- 
element  Reticon  linear  photodiode  array.  The  shift  in 
peak  heights  (owing  to  the  linear  weight)  is  evident, 
indicating  good  qualitative  agreement  with  the  com¬ 
puter-simulation  results. 

This  method  investigates  an  approach  based  on  the 
special  properties  of  Bessel  functions.  If  the  input  / (r ) 
is  correlated  with  Jm(pr)  then  the  Hankel  coefficient 


(a)  (b) 


Fig.  5.  The  oscilloscope  traces  of  the  output  of  the  optical 
processor  performing  the  zeroth-order  Hankel  transform. 
The  traces  correspond  to  |/o(y)|2  versus  y  for  the  two  inputs, 
(a)  For  input  /j  (r)  =  PiJ0(Pir).  (b)  For  input  /2(r)  = 
where  p  >  =  2p\. 


of  order  m,  Fm  (p),  is  obtained  at  the  origin  of  the  cor¬ 
relation  plane. 

Fm(p)  =  f~Kr)Jm[p(r  +  r')]rdr|^o.  (7) 

The  Neumann  addition  theorem  for  Bessel  functions 
states  that 

Jm[p(r  +  r')\=  £  Jm-n(pr')Jn(pr).  (8) 

n*— ® 

Substituting  for  Jm  | p(r  +  r')]  from  Eq.  (8)  into  Eq.  (7) 
we  get 

f(r)Jm[p(r  +  r')]rdr 

=  f  fir)  Z  Jm-n(pr')Jn(pr)rdr 

%/(l 

=  £  Jm.n(pr')Fn(p).  (9) 


Fig.  3.  The  schematic  diagram  of  the  experimental  setup  for 
performing  the  one-dimensional  correlation  between  the 
linearly  weighted  and  coordinate-transformed  input,  fix),  and 
similarly  weighted  and  transformed  Bessel  function  U). 
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Fig.  4.  The  results  of  the  computer  simulation  of  zeroth- 
order  Hankel  transform  of  p\JoU>\r)  and  />.-/, The 
linearly  weighted  and  coordinate  transformed  Hankel  coef¬ 
ficients,  | Fn(y ) |2,  is  plotted  versusy.  Here  *  2/»|. 


Thus,  it  is  seen  that  the  correlation  plane  contains  a  sum 
of  Bessel  functions  of  different  order  weighted  by 
Hankel  transform  coefficients  evaluated  at  p,  of  corre¬ 
sponding  order.  In  principle  it  is  possible  to  extract 
Hankel  transform  coefficients  of  different  order  out  of 
a  single  one-dimensional  correlation  operation.  The 
other  dimension  of  a  two-dimensional  optical  processor 
can  be  used  to  perform  correlations  with  Jm(pr)  with 
different  values  of  p. 
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